home *** CD-ROM | disk | FTP | other *** search
/ SGI Hot Mix 17 / Hot Mix 17.iso / HM17_SGI / research / lib / project_vol.pro < prev    next >
Text File  |  1997-07-08  |  13KB  |  332 lines
  1. ; $Id: project_vol.pro,v 1.8 1997/01/15 03:11:50 ali Exp $
  2. ;
  3. ; Copyright (c) 1992-1997, Research Systems, Inc. All rights reserved.
  4. ;    Unauthorized reproduction prohibited.
  5. ;
  6. ;+
  7. ; NAME:
  8. ;    PROJECT_VOL
  9. ;
  10. ; PURPOSE:
  11. ;       This function returns a two dimensional image that is the
  12. ;       projection of a 3-D volume of data onto a plane (similar to an
  13. ;       X-ray). The returned image is a translucent rendering of the
  14. ;       volume (the highest data values within the volume show up as the
  15. ;       brightest regions in the returned image). Depth queing and
  16. ;       opacity may be used to affect the image. The volume is
  17. ;       projected using a 4x4 matrix, so any type of projection may
  18. ;       be used including perspective. Typically the system viewing
  19. ;       matrix (!P.T) is used as the 4x4 matrix.
  20. ;
  21. ;    Also, see the intrinsic procedure VOXEL_PROJ which performs many
  22. ;    of the same functions as this routine.   VOXEL_PROJ provides the
  23. ;       RGBO rendering method.   VOXEL_PROJ is faster, but it does not
  24. ;       allow for perspective projections.
  25. ;
  26. ;       PROJECT_VOL can combine the contents of the Z-buffer with the
  27. ;       projection when the Z_BUFFER keyword is set.   PROJECT_VOL will not,
  28. ;       however, modify the contents of the Z-buffer.
  29. ;
  30. ; CATEGORY:
  31. ;    Volume Rendering.
  32. ;
  33. ; CALLING SEQUENCE:
  34. ;       Image = PROJECT_VOL(Vol, X_sample, Y_sample, Z_sample)
  35. ;
  36. ; INPUTS:
  37. ;       Vol:        The three dimensional volume of data to project.
  38. ;                   Data type : Any 3-D array except string or structure.
  39. ;       X_sample:   The number of rays to project along the X dimension
  40. ;                   of the image. (The returned image will have the
  41. ;                   dimensions X_sample by Y_sample).
  42. ;                   Data type : Long.
  43. ;       Y_sample:   The number of rays to project along the Y dimension
  44. ;                   of the image.
  45. ;                   Data type : Long.
  46. ;       Z_sample:   The number of samples to take along each ray.
  47. ;                   Higher values for X_sample, Y_sample, and Z_sample
  48. ;                   increase the image resolution as well as execution time.
  49. ;                   Data type : Long.
  50. ;
  51. ; KEYWORD PARAMETERS:
  52. ;       AVG_INTENSITY:
  53. ;                   If set, then the average intensity method of projection
  54. ;                   is used.   The default is a maximum intensity projection.
  55. ;                   This keyword only has effect when NOT using the Z-buffer.
  56. ;       CUBIC:      If set, then the cubic method of interpolation is used.
  57. ;                   The default is a bilinear interpolation.
  58. ;       DEPTH_Q:    Set this keyword to indicate that the image should be
  59. ;                   created using depth queing. The depth queing should
  60. ;                   be a single floating point value (between 0.0 and 1.0).
  61. ;                   This value specifies the brightness of the farthest
  62. ;                   regions of the volume relative to the closest regions
  63. ;                   of the volume. A value of 0.0 will cause the back
  64. ;                   side of the volume to be completely blacked out,
  65. ;                   while a value of 1.0 indicates that the back side
  66. ;                   will show up just as bright as the front side.
  67. ;                   The default is 1.0 (indicating no depth queing).
  68. ;                   Data type : Float.
  69. ;       OPAQUE:     A 3-D array with the same size and dimensions as Vol.
  70. ;                   This array specifies the opacity of each cell in the
  71. ;                   volume. Opaque values of 0 allow all light to
  72. ;                   pass through. Opaque values are cumulative.
  73. ;                   For example, if a ray eminates from a data value of 50,
  74. ;                   and then passes through 10 opaque cells (each with a
  75. ;                   data value of 0 and an opacity value of 5) then that
  76. ;                   ray would be completely blocked out (the cell with the
  77. ;                   data value of 50 would be invisible on the returned
  78. ;                   image).   The default is no opacity.
  79. ;                   Data type : Any 3-D array except string or structure
  80. ;                               (usually the same type as Vol).
  81. ;       TRANS:      A 4x4 floating point array to use as the
  82. ;                   transformation matrix when projecting the volume.
  83. ;                   The default is to use the system viewing matrix (!P.T).
  84. ;                   Data type : Fltarr(4, 4).
  85. ;       XSIZE:      The X size of the image to return.   Congrid is used to
  86. ;                   resize the final image to be XSIZE by YSIZE.   The default
  87. ;                   is the X size of the current window (or the X size of the
  88. ;                   Z-buffer).   If there is no current window then the
  89. ;                   default is X_sample.
  90. ;                   Data type: Int or Long.
  91. ;       YSIZE:      The Y size of the image to return.   Congrid is used to
  92. ;                   resize the final image to be XSIZE by YSIZE.   The default
  93. ;                   is the Y size of the current window (or the Y size of the
  94. ;                   Z-buffer).   If there is no current window then the
  95. ;                   default is Y_sample.
  96. ;                   Data type: Int or Long.
  97. ;       Z_BUFFER:   If set, then the projection is combined with the contents
  98. ;                   of the Z-buffer.   The default is to not use the Z-buffer
  99. ;                   contents.
  100. ;
  101. ; OUTPUTS:
  102. ;       This function returns the projected volume as a two dimensional
  103. ;       array with the same data type as Vol. The dimensions of the
  104. ;       returned array are XSIZE by YSIZE.
  105. ;
  106. ; EXAMPLE:
  107. ;       Use "T3D" to set up a viewing projection and render a volume of
  108. ;       data using "PROJECT_VOL" :
  109. ;
  110. ;       ; Create some data.
  111. ;          vol = RANDOMU(s, 40, 40, 40)
  112. ;          FOR i=0, 10 DO vol = SMOOTH(vol, 3)
  113. ;          vol = BYTSCL(vol[3:37, 3:37, 3:37])
  114. ;          opaque = RANDOMU(s, 40, 40, 40)
  115. ;          FOR i=0, 10 DO opaque = SMOOTH(opaque, 3)
  116. ;          opaque = BYTSCL(opaque[3:37, 3:37, 3:37], TOP=25B)
  117. ;
  118. ;       ; Set up the view.
  119. ;          xmin = 0 & ymin = 0 & zmin = 0
  120. ;          xmax = 34 & ymax = 34 & zmax = 34
  121. ;          !X.S = [-xmin, 1.0] / (xmax - xmin)
  122. ;          !Y.S = [-ymin, 1.0] / (ymax - ymin)
  123. ;          !Z.S = [-zmin, 1.0] / (zmax - zmin)
  124. ;          T3D, /RESET
  125. ;          T3D, TRANSLATE=[-0.5, -0.5, -0.5]
  126. ;          T3D, SCALE=[0.7, 0.7, 0.7]
  127. ;          T3D, ROTATE=[30, -30, 60]
  128. ;          T3D, TRANSLATE=[0.5, 0.5, 0.5]
  129. ;          window, 0, xsize=512, ysize=512
  130. ;
  131. ;       ; Generate and display the image.
  132. ;          img = PROJECT_VOL(vol, 64, 64, 64, DEPTH_Q=0.7, $
  133. ;                OPAQUE=opaque, TRANS=(!P.T))
  134. ;          TVSCL, img
  135. ;
  136. ; MODIFICATION HISTORY:
  137. ;     Written by:    Daniel Carr. Tue Sep  1 17:52:06 MDT 1992
  138. ;
  139. ;       Modified to increase speed.   Also modified
  140. ;       to use current data->normal coordinate conversion.   Added
  141. ;       CUBIC, AVG_INTENSITY, XSIZE, YSIZE, and Z_BUFFER keywords.
  142. ;                       Daniel Carr. Tue Nov 15 16:03:15 MST 1994
  143. ;-
  144.  
  145. FUNCTION Project_Vol, vol, x_sample, y_sample, z_sample, $
  146.          Depth_q=depth_q, Opaque=opaque, Trans=trans, Cubic=cubic, $
  147.          Xsize=xsize, Ysize=ysize, Avg_Intensity=avg_intensity, $
  148.          Z_Buffer=z_buffer
  149.  
  150. ; *** Test inputs.
  151.  
  152. size_vol = Size(vol)
  153. vol_type = size_vol[size_vol[0]+1]
  154.  
  155. x_sample = Long(x_sample[0])
  156. y_sample = Long(y_sample[0])
  157. xy_sample = x_sample * y_sample
  158. z_sample = Long(z_sample[0])
  159. zf_sample = Float(z_sample)
  160. zf_sample_m1 = zf_sample - 1.0
  161. z_sample_m1 = z_sample - 1L
  162. z_max = Float(z_sample - 1L)
  163.  
  164. IF (n_elements(xsize) LE 0L) THEN BEGIN
  165.    IF (!D.Window GE 0L) THEN xsize = !D.X_Size ELSE xsize = x_sample
  166. ENDIF
  167.  
  168. IF (n_elements(ysize) LE 0L) THEN BEGIN
  169.    IF (!D.Window GE 0L) THEN ysize = !D.Y_Size ELSE ysize = y_sample
  170. ENDIF
  171.  
  172. IF (N_ELEMENTS(depth_q) GT 0) THEN depth_q = (Float(depth_q[0]) > 0.0) < 1.0 $
  173. ELSE depth_q = 1.0
  174. depth_q = 1.0 - depth_q
  175.  
  176. block_out = 0B
  177. IF (N_Elements(opaque) GT 0L) THEN BEGIN
  178.    IF (N_Elements(opaque) EQ N_Elements(vol)) THEN BEGIN
  179.       opaque = Reform(opaque, size_vol[1], size_vol[2], size_vol[3])
  180.       block_out = 1B
  181.    ENDIF ELSE BEGIN
  182.       Print, 'Opaque array must be the same size as volume array.'
  183.       RETURN, Bytarr(xsize, ysize)
  184.    ENDELSE
  185. ENDIF
  186.  
  187. IF (N_Elements(trans) GT 0) THEN BEGIN
  188.    IF (N_Elements(trans) NE 16) THEN BEGIN
  189.       Print, 'Incorrect number of elements in Trans.'
  190.       RETURN, Bytarr(xsize, ysize)
  191.    ENDIF ELSE BEGIN
  192.       trans = Float(Reform(trans, 4, 4))
  193.    ENDELSE
  194. ENDIF ELSE BEGIN
  195.    trans = !P.T
  196. ENDELSE
  197. trans = Invert(trans, status)
  198. IF (status NE 0) THEN BEGIN
  199.    Print, 'Unable to invert transformation matrix.'
  200.    RETURN, Bytarr(xsize, ysize)
  201. ENDIF
  202.  
  203. ; *** Set up the required variables.
  204.  
  205. x_ind = (((Float(size_vol[1]) * Findgen(x_sample) / Float(x_sample - 1L)) * $
  206.    !X.S[1]) + !X.S[0]) # Replicate(1.0, y_sample)
  207. y_ind = Replicate(1.0, x_sample) # $
  208.         (((Float(size_vol[2]) * Findgen(y_sample) / Float(y_sample - 1L)) * $
  209.    !Y.S[1]) + !Y.S[0])
  210.  
  211. max_vol = Float(Max(vol))
  212.  
  213. IF Keyword_Set(z_buffer) THEN BEGIN ; *** Use Z-Buffer info.
  214.    save_win = !D.Window
  215.    save_name = !D.Name
  216.    Set_Plot, 'Z'
  217.    img = Congrid(Tvrd(), x_sample, y_sample, /Interp, /Minus_One)
  218.    min_img = 0B
  219.    depth_fac = 0B
  220.    depth = Congrid(Tvrd(Channel=1, /Words), x_sample, y_sample, $
  221.       /Interp, /Minus_one)
  222.    xsize = !D.X_size
  223.    ysize = !D.Y_size
  224.    Set_Plot, save_name
  225.    Wset, save_win
  226.  
  227.    ; *** Do the projection.
  228.  
  229.    FOR k=0L, z_sample_m1 DO BEGIN
  230.       kf = Float(k)
  231.       z_pos = (kf * !Z.S[1]) + !Z.S[0]
  232.  
  233.       index = [[x_ind[*]], [y_ind[*]], $
  234.          [replicate(z_pos, xy_sample)], [replicate(1.0, xy_sample)]] # trans
  235.  
  236.       indx = Float(size_vol[1]) * index[*, 0] / index[*, 3]
  237.       indy = Float(size_vol[2]) * index[*, 1] / index[*, 3]
  238.       indz = Float(size_vol[3]) * index[*, 2] / index[*, 3]
  239.  
  240.       z_pos = 2.0 * ((((z_pos - 0.5) > (-0.5)) < 0.5) * 32765.0)
  241.       z_ind = Where(z_pos GT depth, count)
  242.  
  243.       IF (count EQ (xsize * ysize)) THEN BEGIN ; *** Nothing in the way.
  244.          IF (block_out) THEN $
  245.             img = img - $
  246.                (Interpolate(opaque, indx, indy, indz, Missing=min_img, $
  247.                             cubic=Keyword_Set(cubic)) < img)
  248.  
  249.          depth_fac[0] = ((1.0 - (kf / z_max)) * depth_q * max_vol)
  250.          img = img > ((Interpolate(vol, indx, indy, indz, Missing=min_img, $
  251.                                    cubic=Keyword_Set(cubic)) > $
  252.             depth_fac[0]) - depth_fac[0])
  253.       ENDIF ELSE BEGIN ; *** Something in the way.
  254.          IF (count GE 1L) THEN BEGIN ; *** Something new in front.
  255.             IF (block_out) THEN BEGIN
  256.                temp_img = $
  257.                   Interpolate(opaque, indx, indy, indz, Missing=min_img, $
  258.                               cubic=Keyword_Set(cubic)) < img
  259.                img[z_ind] = img[z_ind] - temp_img[z_ind]
  260.             ENDIF
  261.  
  262.             depth_fac[0] = ((1.0 - (kf / z_max)) * depth_q * max_vol)
  263.             temp_img = $
  264.                (Interpolate(vol, indx, indy, indz, Missing=min_img, $
  265.                             cubic=Keyword_Set(cubic)) > $
  266.                depth_fac[0]) - depth_fac[0]
  267.             img[z_ind] = img[z_ind] > temp_img[z_ind]
  268.          ENDIF ; *** Nothing new in front.
  269.       ENDELSE
  270.    ENDFOR
  271. ENDIF ELSE BEGIN ; *** Don't use Z-Buffer.
  272.    IF Keyword_Set(avg_intensity) THEN img = Fltarr(x_sample, y_sample) $
  273.    ELSE BEGIN
  274.       CASE vol_type OF
  275.          2: img = Intarr(x_sample, y_sample)
  276.          3: img = Lonarr(x_sample, y_sample)
  277.          4: img = Dblarr(x_sample, y_sample)
  278.          5: img = Fltarr(x_sample, y_sample)
  279.          6: img = Complexarr(x_sample, y_sample)
  280.       ELSE: img = Bytarr(x_sample, y_sample)
  281.       ENDCASE
  282.    ENDELSE
  283.  
  284.    min_img = img[0]
  285.    depth_fac = img[0]
  286.  
  287.    ; *** Do the projection.
  288.  
  289.    FOR k=0L, z_sample_m1 DO BEGIN
  290.       kf = Float(k)
  291.       z_pos = (kf * !Z.S[1]) + !Z.S[0]
  292.  
  293.       index = [[x_ind[*]], [y_ind[*]], $
  294.          [replicate(z_pos, xy_sample)], [replicate(1.0, xy_sample)]] # trans
  295.  
  296.       indx = Float(size_vol[1]) * index[*, 0] / index[*, 3]
  297.       indy = Float(size_vol[2]) * index[*, 1] / index[*, 3]
  298.       indz = Float(size_vol[3]) * index[*, 2] / index[*, 3]
  299.  
  300.       IF (block_out) THEN $
  301.          img = img - $
  302.          (Interpolate(opaque, indx, indy, indz, Missing=min_img, $
  303.                       cubic=Keyword_Set(cubic)) < img)
  304.  
  305.       depth_fac[0] = ((1.0 - (kf / z_max)) * depth_q * max_vol)
  306.       IF Keyword_Set(avg_intensity) THEN $
  307.          img = img + (((Interpolate(vol, indx, indy, indz, Missing=min_img, $
  308.                         cubic=Keyword_Set(cubic)) > $
  309.             depth_fac[0]) - depth_fac[0]) / zf_sample) $
  310.       ELSE $
  311.          img = img > ((Interpolate(vol, indx, indy, indz, Missing=min_img, $
  312.                        cubic=Keyword_Set(cubic)) > $
  313.             depth_fac[0]) - depth_fac[0])
  314.    ENDFOR
  315.    IF Keyword_Set(avg_intensity) THEN BEGIN
  316.       CASE vol_type OF
  317.          2: img = Fix(img)
  318.          3: img = Long(img)
  319.          4: img = Double(img)
  320.          5: img = Float(img)
  321.          6: img = Complex(img)
  322.       ELSE: img = Byte(img)
  323.       ENDCASE
  324.    ENDIF
  325. ENDELSE
  326.  
  327. IF ((xsize NE x_sample) OR (ysize NE y_sample)) THEN $
  328.   img = Congrid(img, xsize, ysize, /Interp, /Minus_One)
  329.  
  330. RETURN, img
  331. END
  332.